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AN  APPLICATION  OF  THE  METHOD  OF  LINES 
TO  THE  TRANSONIC  AIRFOIL  PROBLEM 

I .  Introduction 

Background 

The  demands  for  faster  and  more  efficient  subsonic 
aircraft  have  forced  designers  to  reconsider  the  design 
problems  associated  with  transonic  flow.  Until  recently, 
the  complications  inherent  in  this  flight  realm  could  be 
avoided  by  designing  the  aircraft  for  its  optimum  perfor¬ 
mance  either  below  or  above  the  transonic  regime.  With  this 
approach  the  surrounding  flow  field  could  be  considered  as 
being  either  entirely  subsonic  or  completely  supersonic; 
however,  these  simplified  approaches  are  impossible  in  the 
design  of  the  aircraft  presently  being  asked  for.  Cargo 
and  passenger  planes  that  can  fly  at  very  high  subsonic 
speeds  and  interceptors  required  to  fly  at  these  same  speeds 
to  most  effectively  deliver  their  weapons  will  want  their 
optimum  periormance  in  the  transonic  portion  of  the  velocity 
spectrum,  'lii^refore,  solutions  to  the  long  standing  transonic 
flow  problem  are  now  a  necessity  and  are  of  prime  importance 
in  this  the  initial  stage  of  research  and  development. 

The  first  and  perhaps  most  difficult  hurdle  that  must 
be  crossed  in  transonic  design  problems  is  the  development 
of  a  solution  technique  for  some  mathematical  model  of  the 
flow  field  about  a  2-dimensional  airfoil  in  transonic  flight. 
There  are  several  reasons  why  acrodynamicists  have  found  this 
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to  be  difficult  and  have  avoided  it  whenever  possible;  all 
of  the  reasons  are  related  to  what  is  physically  occuring 
to  the  flow  in  the  region  of  the  airfoil.  By  definition,  the 
transonic  airfoil  problem  is  concerned  with  determining  the 
pressure  distribution  over  an  airfoil  when  the  flow  field 
around  it  is  a  mixture  of  subsonic  (M  <  1)  and  supersonic 
(M  >  1)  flow  regions.  When  these  regions  coexist  in  a  field 
it  is  referred  to  as  a  supercritical  flow  while  a  subcritical 
flow  is  defined  to  be  a  flow  in  which  the  velocity  is  sub¬ 
sonic  at  every  point  and  a  critical  flow  refers  to  flow  when 
M  s  1  at  only  one  point  in  the  field.  There  are  two  important 
mathematical  complications  that  arise  when  the  flow  is 
supercritical.  First,  if  the  third  independent  variable, 
namely  time,  is  not  retained  in  the  differential  equations, 
then  the  basic  form  of  the  equations  changes  from  elliptic 
in  subsonic  regions  to  hyperbolic  in  supersonic  regions. 
Second,  for  both  the  steady  and  unsteady  problem,  as  the 
freestream  Mach  number  approaches  unity,  nonlinear  terms 
must  be  retained  in  the  small  perturbation  equations  often 
used  as  a  mathematical  model.  In  addition,  a  third  major 
mathematical  difficulty  appears  in  attempting  to  model  the 
effects  of  the  shock  wave  that  is  normally  attached  to  the 
airfoil  surface.  To  account  foi  this  shock  either  a  compli¬ 
cated  system  of  conservation  equations  in  divergence  form 
must  be  used  or  additional  relationships  (Rankine-Hu" oniot 
equations)  must  be  included  and  thus  amplify  the  complexity 
of  the  problem.  With  all  these  added  difficulties  induced 
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by  the  physical  nature  of  the  transonic  flow  field,  the 
mathematical  problem  has  defied  exact  solution  except  for 
extremely  limited  cases.  Fortunately,  though,  the  develop¬ 
ment  and  availability  of  high  speed  computers  has  produced 
approximate  solutions  by  numerical  techniques  which  in  some 
cases  are  in  better  agreement  with  flight  data  than  wind  tunnel 
data. 

Statement  of  the  Problem 

The  objective  of  this  work  is  to  develop  and  evaluate 
a  particular  numerical  solution  technique  for  the  two- 
dimensional,  inviscid  transonic  airfoil  problem.  The  method 
is  required  to  serve  as  an  efficient  and  workable  design  tool 
for  the  evaluation  of  ar  airfoil  in  its  transonic  regime; 
therefore,  it  must  yield  acceptable  results  with  minimum 
computer  time  and  storage  space.  As  presented  here  the 
method  is  applied  only  to  thin,  symmetrical  airfoils  with 
sharp  leading  edges.  This  permits  the  present  study  to 
concentrate  on  developing  and  testing  the  solution  technique 
and  then,  if  warranted,  future  efforts  could  be  concentrated 
on  extending  the  technique  to  problems  of  more  extensive 
scope.  Nothing  will  be  included  to  explicitly  account  for 
the  presence  of  shock  waves  in  the  flow,  and  the  associated 
viscous  interaction  with  the  boundary  layer  will  not  be 
considered.  This  last  assumption,  though,  is  made  simply 
because  it  is  required  to  reduce  the  mathematical  problem 
to  the  scope  of  this  presentation;  the  intent  is  not  to 
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imply  that  inviscid  flow  is  always  a  good  assumption  for  this 
velocity  regime.  Pearcy  (Ref  -)  has  shown  that  the  conpli- 
cated  interaction  between  the  shock  wave  and  the  boundary 
layer  will  result  in  flow  separation  and  disturbances  not 
present  in  flow  fields  outside  the  transonic  regime. 

Therefore,  it  is  agreed  that  these  viscous  effects  could 
indeed  be  a  dominant  feature  for  many  practical  problens; 
however,  at  this  time  there  is  no  practical  way  to  approach 
the  full  viscous  problem  within  the  scope  proposed  here. 

It  is  hoped  that  the  solution  technique  developed  in  this 
presentation  nay  be  extendable  to  the  viscous  problen  but 
the  method  as  presented  here  will  only  be  applicable  to 
transonic  airfoil  problems  where  the  viscous  effects  do  not 
drastically  alter  the  flow  field. 

Other  Solutions  Currently  Available 

The  approach  to  this  problem  was  selected  after 
examining  existing  solution  methods.  As  implied  before, 
exact  solution  techniques  for  the  inviscid  transonic  airfoil 
problem  are  available;  however,  they  contain  well-nigh 
insurmountable  difficulties  except  for  extremely  limited 
cases  as  shown  by  Ferrari  and  Tricomi  (Ref  2:562).  Because 
of  the  lack  of  exact  solutions  many  numerical  methods  have 
been  devised  to  produce  approximate  solutions.  These,  in 
general,  may  be  classified  into  one  of  three  broad  categories. 
In  the  first  category,  the  governing  equations  are  trans¬ 
formed  into  integral  relations  and  solved  by  iterative 
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procedures.  This  net hod  h2S  been  employed  by  Osvatitsch 
(Ref  2:563-574)  and  more  recently  by  Spreiter  and  Alksne 
(Ref  10)  in  solving  the  transonic  problem  for  thin,  sharp¬ 
nosed  airfoils.  Although  this  technique  could  satisfy  the 
computation  time  objective,  at  the  moment  it  is  not  extend¬ 
able  to  arbitrary  airfoil  shapes.  The  second  category  of 
solutions  transforms  the  differential  equations  into  alge¬ 
braic  expressions  by  an  explicit  or  implicit  finite-difference 
scheme.  Magnus  and  Yoshihara  (Ref  7)  have  been  particularly 
successful  with  this  technique  and  can  now  produce  extremely 
good  results  for  several  airfoil  shapes.  The  difficulty 
associated  with  these  methods  is  the  long  computer  times 
(on  the  order  of  1  hour/input  for  the  Magnus  2nd  Yoshihara 
method)  and  storage  space  required.  This  nakes  the  method 
unattractive  as  a  design  tool  and  infeasible  to  extend  to 
the  more  involved  problems.  The  third  category  of  solution 
techniques  reduces  the  system  of  governing  partial  differen¬ 
tial  equations  to  a  set  of  ordinary  differential  equations 
which  may  be  solved  by  one  of  several  available  solution 
procedures.  Dorodnitsyn  (Ref  1)  originally  suggested  this 
approach  to  the  transonic  problem,  and  the  recent  successful 
application  of  these  methods  by  Tai  (Ref  11)  and  Melnik  2nd 
Ives  (Ref  8)  prompted  the  decision  to  use  a  version  of  this 
technique  for  the  development  of  the  proposed  method.  By 
taking  a  very  simplified  method  of  weighted  residuals 
approach,  Tai  has  shown  ’’hat  this  method  is  not  only  appli¬ 
cable  to  the  transonic  airfoil  problem,  but  also  that  it  can 
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produce  excellent  results  Kith  very  little  computer  time. 
Meanwhile,  Melnik  and  Ives  have  used  a  similar  approach  and 
have  shown  that  weighting  functions  of  the  subdomain  type  are 
capable  of  handling  any  arbitrary  airfoil  shape  if  the  flow 
field  is  of  finite  domain. 

Approach  to  the  Problem 

Small  perturbation  theory  2nd  irrotationality  will  be 
used  to  provide  the  set  of  governing  differential  equations. 
The  two  boundary  conditions  these  equations  will  be  subject 
to  are:  no  flow  normal  to  the  airfoil  surface  and  undis¬ 
turbed  flow  at  infinity.  Prior  to  developing  the  solution 
technique,  the  fora  of  the  governing  equations  will  be 
changed  by  nondiaensionalizing  the  variables  and  making  a 
coordinate  transformation  to  reduce  the  domain  to  a  finite 
region.  This  is  done  to  sinplify  application  of  the  numerical 
technique  and  to  eliminate  the  difficulty  associated  with 
applying  the  boundary  condition  at  infinity. 

The  particular  solution  technique  to  be  applied  is  a 
form  of  the  nethod  of  lines  to  reduce  the  set  of  governing 
equations  to  a  set  of  ordinary  first-order  differential 
equations.  This  new  set  is  then  solved  by  the  classical 
Runge-Kutta  fourth-order  technique  for  the  perturbation 
velocities  throughout  the  field.  Finally,  the  solution 
technique  is  evaluated  by  comparing  the  analytical  results 
with  published  experimental  data  for  symmetrical  airfoils 
at  zero  incidence  angle. 
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II .  Formulation  of  the  Solution 
Technique 

In  developing  an  approach  to  the  transonic  airfoil 
problem  three  very  basic  but  crucial  decisions  must  be  made. 
First,  a  set  of  governing  equations  must  be  selected  which 
completely  define  the  problem  within  some  prescribed  scope. 
Next,  modifications  of  the  set,  such  as  coordinate  or  variable 
transformations  must  be  considered  if  they  will  enhance  the 
solution  or  simplify  the  mechanics  of  obtaining  the  solution. 
Finally,  from  the  variety  of  available  numerical  schemes, 
one  must  be  selected  to  solve  the  resulting  system  of 
equations.  The  product  of  these  three  decisions  defines  a 
specific  solution  technique;  therefore,  the  alternatives  to 
each  of  these  choices  will  be  discussed  in  explaining  the 
rationale  for  the  approach  used  in  this  presentation. 

Picking  a  set  of  governing  equations  for  a  solution 
method  requires  a  choice  to  be  cade  based  on  the  exactness 
demanded  of  the  solution  versus  the  amount  of  mathematical 
complexity  acceptable  in  the  method.  For  inviscid  flow 
there  are  three  different  sets  of  relationships  which  are 
commonly  used  to  define  the  transonic  airfoil  problem  and 
the  complexity  of  each  of  these  sets  is  directly  related  to 
how  the  shock  waves  are  mathematically  accounted  for.  The 
first  set  consists  of  the  iscntropic  relationship 
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and  tee  unsteady  form  of  the  conservation  equations  for  mass 
and  m omen turn 
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These  are  used  throughout  the  flow  field  except  across  the 
shock  where  the  additional  Rankine-Hugoniot  relationships 
are  required-  Grossman  and  Horetti  (Ref  4)  used  this  choice 
in  their  solution  technique  but  this  decision  resulted  in 
adding  appreciable  comp 15 cations  when  the  method  was  applied 
to  transonic  airfoils.  The  second  alternative  is  to  use  the 
steady  state  fora  of  the  same  set  of  equations  but  written 
in  divergence  fora  because  in  this  form  the  system  is  then 
applicable  even  across  the  shock  wave.  This  set  still 
consists  of  three  nonlinear  partial  differential  equations 
together  with  an  algebraic  relationship  and  therefore 
recains  a  very  conplex  sys ten  to  work  with.  Although  Tai 
(Ref  11)  developed  a  solution  technique  based  on  fhe  sub- 
donain  octhod  using  this  set  of  governing  equations,  a 
sinpler  systen  is  desirable  for  this  presentation.  A  third 
possible  choice  exists  when  it  is  assuned  that  a  good  approx- 
ination  to  ..he  flow  field  can  be  obtained  without  nathenati- 
cally  accounting  for  the  shock  waves.  Sprcitcr  and  Alksne 
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(Ref  10)  2nd  several  others  save  used  this  2Sstmption  while 
working  with  integral  solutions  for  slender,  sharp-nosed 
bodies  in  the  transonic  regime.  This  choice  seems  particularly 
attractive  for  this  presentation  since  the  scope,  as  stated 
in  the  introduction,  has  already  excluded  the  interaction 
between  the  shock  and  the  boundary  layer  and  becaese,  ia 
general,  no  a  priori  conclusion  can  be  drawn  concerning  the 
magnitude  of  the  change  in  the  flow  due  exclusively  to  the 
imriscid  shock  effects.  Further,  since  this  is  a  preliminary 
study  of  a  new  solution  technique  this  assumption  permits  the 
evaluation  of  the  technique  with  the  minimum  mathematical 
complexity.  Then,  once  the  technique  is  formulated  using 
the  simpler  set,  the  same  principles  could  be  applied  to  the 
more  exact  equations  if  the  solution  required  it. 

Kith  this  decision,  small  perturbation  theory  and  the 
irrotationality  condition  can  be  used  to  define  a  set  of 
governing  equations-  Liepmann  and  Roshko  (Ref  6:202-2OS) 
develop  the  theory  and  give  the  form  of  the  small  perturba¬ 
tion  relationship  which  is  valid  throughout  the  transonic 
regine: 


(1 


«  2^  3u  .  3v  _  „  2  Y  *  1 
3x  37  “  ~  ~ - 


CO 


(1) 


This  equation  when  nondinensionali zed  and  rearranged 
(Appendix  A)  nay  be  written 


1(1  -  M„2)«  -  \  Mm2(Y  -  Du2)  ♦  !y  =  0  (2) 
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Traditionally ,  the  coalitiearity  of  Eq  (2)  and  t he  change  is 
sign  of  the  first  term  for  some  positive  value  of  vs  as 
approaches  one  has  made  this  equation  difficult  to  solve, 
la  an  attempt  to  relocate  these  difficulties  into  a  form 
which  might  he  more  easily  dealt  with  numerically  the 
following  definition  is  made 


2  =  (i  -  -  4  M»2(?  +  *)*2 

2 


Therefore  the  fall  set  of  governing  equations  needed  to 
completely  define  the  problem  consists  of  Eq'  (3).  the  modi¬ 
fied  small  perturbation  relationship. 


3  C  3v 

♦  — —  =  o 
3x  3y 


and  the  condition  of  irrotationali ty 


~  | H-  =  0  C5) 

3x  3y 

Further,  these  partial  differential  equations  are  subject  to 
boundary  conditions  on  the  body  surface  and  at  infinity.  If 
the  airfoil  surface  is  defined  by  the  relation  v  =  fix'), 
then  the  condition  at  the  body  surface  of  no  flow  normal  to 
the  surface  (Ref  6:208)  nay  be  written 


u4-U^  *  37  =  *’  <6> 

Because  snail  disturbance  theory  is  being  used  for  the 
governing  equations  the  only  nathcnatical  boundary  condition 
that  is  specified  at  infinity  is  that  u  and  v  are 
both  finite.  However,  at  infinity  in  the 


10 


6AM/AE/72-7 


physics!  probleat  o  2nd  v  both  sre  zero  and  this  more 
restrictive  condition  is  generally  assumed  whenever  there 
are  so  regions  of  supersonic  flow  extending  to  the  infinity 
boundary-  Therefore  for  this  solution  technique  the  boundary 
condition  which  will  be  imposed  at  infinity  is 


c  =  v  =  £2  =  0 


(7) 


Based  on  the  difficulty  experienced  by  other  authors 
in  applying  the  boundary  conditions,  it  was  decided  that  a 
coordinate  transformation  should  be  used  to  avoid  assumptions 
as  to  where  the  boundary  conditions  should  be  applied.  Sew 
independent  variables  were  defined  as 


^  = 


|x|  *  L 


(8) 


=  L-i-lizl 
h\  *  l 


(9) 


where  L  is  soae  positive  constant.  As  can  be  seen  from 
Fig.  1  and  2  on  the  following  page,  this  transformation 
accomplishes  two  important  objectives.  First,  the  domain 
of  the  independent  variables  is  reduced  to  a  finite  region; 
therefore,  the  boundary  conditions  at  infinity  in  the  x-y 
plane  are  applied  at  unity  in  the  £-n  plane.  Second,  the 
thickness  of  any  airfoil  shape  is  reduced  to  zero  in  the  new 
coordinate  systen  such  that  the  boundary  condition  on  the 
airfoil  surface  is  applied  at  n  =  0. 


1 1 
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Fig.  i.  Original  Reference 
Frame. 


Fig-  2. 


Trass rormed  Rererence 
Frame. 


This  transformation  is  thoroughly  discussed  and  developed 
in  Appendix  A  where  the  small  perturbation  relation  is  shown 
to  be 


32 


Lf 1  (l  -  n) 


30 


L{ 1  -  n) 


H  (L  ♦  f)d-!Ci)2  3n  (L  +  f)(i-K})2  3p 


—  (10) 


and  the  irrotationality  condition  becones 

3v  =  Lf 1  ( 1  -  n) _  3v  +  L(1  -  n)2  3u 

35  (L  +  f)(i-icj)2  (L  *  f)(i-|5|)2  an 

The  algebraic  relationship  renains  unchanged,  therefore, 

Eqs  (3),  (10),  and  (11)  represent  the  full  set  of  transforned 
governing  equations.  The  boundary  conditions  which  these 


equations  arc  subject  to  are 
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u  =  v  =  2  *  0  (12) 

at  C  or  t)  equal  *  1  and 

t  =  r*(l  +  ts) 

at  tiae  body  surface  (si  =  0).  This  last  condition  is  farther 
simplified  by  recalling  that  hr  definition  a  is  much  less 
than  one,  therefore  the  boundary  condition  can  be  approxi¬ 
mated  hr 

r  =  f‘  (13) 

at  ii  =  0. 

The  most  difficult  decision  that  had  to  be  made  was  the 
choice  of  a  numerical  method  to  solve  the  resulting  set  of 
equations.  Because  of  the  objective  to  reduce  computer  tine 
and  storage  space  to  a  minimus  and  due  to  the  successful 
applications  by  Melnik  and  Ives  (Ref  8)  and  Tax,  it  was 
decided  very  early  that  a  form  of  the  Bethod  of  weighted 
residuals  should  be  used.  Hithin  this  approach  several 
alternatives  were  tried  and  a  solution  technique  using  the 
subdomain  Bethod  was  extensively  developed.  Unfortunately, 
the  lack  of  any  rational  procedure  for  initializing  the 
variables  caused  the  nethod  to  be  rejected.  Thus,  although 
a  detailed  discussion  of  this  technique  is  included  as 
Appendix  B,  this  presentation  uses  the  easier  to  apply 
nethod  of  lines. 
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The  purpose  of  the  method  of  lines  is  to  reduce  a 
system  of  partial  differential  equations  to  a  set  of  first- 
order  ordinary  differential  equations.  The  method  of  lines 
presumes  that  the  derivatives  in  the  n-direction  can  be 
represented  by  a  finite-difference  approximation  (details 
are  included  in  Appendix  A).  To  use  this  concept,  the 
tl-domain  is  divided  by  ( !i  -  1)  parallel  lines,  equally 
spaced  between  n  =  0  and  n  =  1 . 


%  =  |  ;  1  =  i,  2, 


...,  (N-l) 


Representing  the  derivatives  with  respect  to  n  by 
central  differences,  the  governing  equations  (3),  (10),  and 
(11)  for  the  ith  line  may  be  written: 


=  (1  -  H„2)ui  -  I  M£02(y  *  l)Ui: 


(14) 


dfi. 


f*(N  -  i ) L 


d£  2  (L  ♦  f)(l  -  U|) 


+  l  '  Oi-x)  ' 


L(N  -  i) 


2K(L  +  f)  (1  -  Kl) 


(vi+l  -  vi-l>  (15) 


dv4 


f ' (N  -  i)  L 


2  (L  +  f)(l  -  Ul) 


(vi  +  l  "  Vi-P 


L(N  -  i)2 


, ' rv”  ^ui+i  -  ui-i)  (16) 

2N(L  +  f) (1  -  \i\)2 
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When  this  system  is  collected  for  all  (N  -  1)  lines  the 
resulting  set  of  (3N  -  3}  equations  involve  (3N  +  3)  unknowns 
u . ,  v. ,  and  12^ ;  j  =  0,  1,  2,  n •  Three  of  the  six 

additional  relationships  which  are  required  to  uniquely 
define  the  problem  are  obtained  from  the  boundary  condition 

uN  =  vN  =  =  0  (17) 

at  n  =  1.  Two  more  relations  are  obtained  by  writing 
governing  equations  (3)  and  (10)  along  the  body  surface 

£0  =  (1  ~  M„2)u0  '  \  +  Du02  (13) 


dfiQ 

d£ 


LN 


2(1  ♦  f)(l  -  Ul) 


2  [f'(-3a0  *  412! 


122)  - 


( "3v0  +  4vj  -  v2)] 


(19) 


where  the  n-deri vatives  are  approximated  by  forward- 
differences.  Finally,  the  last  required  equation  comes 
from  the  boundary  condition  on  the  surface,  Eq  (13),  which 
is  rewritten  as, 


vn  =  f’(x) 


f  • 


Kl 

-  Ul 


/ 


(20) 


Therefore,  the  complete  problem  now  consists  of  2N  first- 
order  ordinary  differential  equations  which  must  be  solved 
simultaneously  with  (N  +  3)  algebraic  relations. 


This  set  can  be  conveniently  solved  if  initial  values 
can  be  determined  for  starting  the  solution  of  the  ordinary 
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differential  equations  and  if  a  unique  value  for  u  can  be 
found  for  a  given  ft  from  Eq  (3)  .  If  Eq  (3)  is  solved  for 
u  by  the  quadratic  formula 


(1  -  Mw2) 

m„2(y  ♦  1) 


1  t 


2Mco2(Y  ♦  1) 
(1  -  M „2)2 


a 


(22) 


then,  a  choice  of  signs  must  be  made  to  uniquely  define  u. 

A  typical  plot  of  u  versus  ft  from  Eq  (3)  for  a  less  than 

one  is  shown  below  in  Fig.  3  where  uc  refers  to  the  value  of 
3  ft 

u  at  which  =  0;  namely. 


(1  -  Mw2) 

=  — 2 - 

M.  (Y  *  1) 


(23) 


and  ft£  is  the  corresponding  maximum  value  of  ft.  Comparing 
Eq  (22)  and  Fig.  1  it  can  be  concluded  that  the  proper  sign 
in.  Eq  (22)  is  negative  if  u  <  uc  and  positive  if  u  >  uc. 


Fig.  3.  Typical  Relation  of  u  to  ft  as  found  from 
Eq  (22)  Including  the  Critical  Values  of  u. 
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Frosi  experimental  data  presented  by  Knechtei  (Kef  S)  it  can 

be  observed  that  for  snbcritical  flows  a  is  always  less  than 

a  while  for  supercritical  flows  u  exceeds  this  value  over 
c 

score  portion  of  the  airfoil.  Typical  examples  of  the 

a-distribation  on  the  surface  of  a  subcritical  (H^  =  0.806} 

and  a  supercritical  =  0.909)  6%  circular  arc  airfoil  are 

given  in  Figs.  4  and  5  below.  Also  shown  in  these  figures 

are  the  values  of  uc  for  each  case  and  the  corresponding 

distribution  of  2  required  when  the  proper  sign  is  used  in 
0  *  ” 

Zq  (22) . 


Mg.  4.  Subcritical 
(Mro  =  .806)  distribution 
of  u,  corresponding 
fi-distribution  from  Eq  (22) 
and  critical  value  of  u  for 
a  6%  circular  arc  airfoil. 


Fig.  5.  Supercritical 
(Moo  =  .909)  distribution 
of  u,  corresponding 
2-distribution  from  Eq  (22) 
and  critical  value  of  u  fo* 
a  6%  circular  arc  airfoil. 
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tShen  the  set  of  governing  eqazt ions  is  solved  by  a 
marching  method  such  2s  Essnge-Eatta,  the  proper  sign  in 
Eq  (22}  can  be  predicted  by  the  value  of  2=  and  2- *  at  the 
points  where  the  value  of  Oj  is  found  to  equal  ec.  There¬ 
fore,  if  the  proper  initial  conditions  of  v^  and  2^  can  be 
determined  £q  (22)  is  solved  for  the  corresponding  a-. 

Then  the  Bonge-Kotta  fourth-order  solution  technique  can  be 
applied  to  solve  the  2K  ordinary  differential  conations  for 
the  value  of  v^  and  some  &£  downstream  and  thereby 
setting  up  the  marching  technique  in  the  ^-direction.  From 
this  nethod  the  values  of  u  and  v  are  £ou«e  at  every  point 
in  a  i?  x  hxi  grid  which  covers  the  entire  flow  field.  With 
these  values  three  important  parameters  of  the  field  are 
found:  the  pressure  distribution  over  the  airfoil  surface 

Cp  =  -  2u0  (24) 

the  ratio  of  local  to  free  strean  velocity  at  each  grid  point 

=  /(ui  +  l)2  ♦  vi2  (25) 

i 

and  the  Mach  number  distribution  throughout  the  f  .eld 

1 

Mi  =  ^MTO[1  -  \  ^  -  ^^Ui  +  Ui2  +  Vi2))  2  (26) 

The  solution  technique  can  then  be  evaluated  by  a  comparison 
with  experimental  data. 
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III .  Bespits 

Tee  resalts  aadl  conclusions  obtained  taroogfe  tbe  coarse 
of  this  stedy  will  be  girea  ia  ccroaologi cal  order  so  that 
the  subject  matter  C2n  be  divided  up  late  fear  main  topics; 
sabcritical  airfoil  steadies,  sapercritical  airfoil  efforts, 
attempts  to  generalize  the  initial  conditions,  and  applica¬ 
tions  of  the  method  to  snbsoaic  and  supersonic  wavy  vails. 

Sabcritical  Airfoi 1 

It  W2S  assumed  In  the  formulation  of  the  solution 
technique  .that  only  symmetrical  airfoils  at  zero  incidence 
angle  would  be  considered.  Therefore,  because  of  avaiiabi lity 
of  experimental  data  accumulated  by  Enechtel  (Ref  5)  the  6% 
circular  arc  airfoil  W2S  selected  2S  the  reference  airfoil 
and  the  suberitical  case  of  M0  =  0.S06  was  chosen  for  the 
first  test  of  the  proposed  solution  technique. 

The  difficulties  in  generating  initial  conditions  for 

upstream  of  the  airfoil  for  this  problem  will  be  discussed 

in  detail  later  in  this  section  but  a  method  was  constructed 

at  this  point  for  determining  the  initial  conditions  by 

starting  the  solution  at  a  point  on  the  airfoil  (0  <  x  £  1). 

Based  on  the  form  of  the  experimental  Cp  distribution  and 

the  curve  of  the  boundary  condition  Vg  =  f’  (as  shown  in 

Fig.  f'  on  the  next  page)  it  was  assumed  that  the  distribution 

of  any  u^(£)  (recalling  that  the  subscripts  on  u,  v,  or  SI 

refer  to  their  value  on  the  line  rj  =  —  )  is  symmetrical  with 

N 

respect  to  the  line  x  =  1/2  while  the  distribution  of  any 
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Fig-  6.  : rp i c s i  Sistribution  of  ug  2nd  Vg  Over 

2  Circular  Arc  Airfoil - 

v-(C)  is  symmetrical  with  respect  to  the  point  x  =  1/2, 
y  =  f(l/2).  With  these  assumptions  the  initial  values  of 
and  V-  can  be  guessed  at  a  point  ahead  of  x  =  1/2,  then 
fcnoving  what  the  values  of  u^  and  should  be  at  some  point 
downstreaa,  these  initial  guesses  are  refined  by  an  iteration 
process.  In  application,  froa  the  experimental  data  the 
point  at  which  Cp  =  0  was  selected  as  the  initial  value  of 
x,  the  values  of  U;  at  this  point  vere  assuned  to  be  zero, 
and  the  values  of  were  guessed.  Since  the  assuned 
symmetry  required  all  to  be  zero  at  x  =  1/2,  these  guesses 
for  the  initial  conditions  were  refined  by  iterating  between 
the  initial  x  and  x  =  1/2  until  all  the  changed  signs 
within  a  region  x  =  0.5  ±  0.002.  When  this  condition  was 
met  it  was  found  that  the  entire  flow  field  (0  <  x  <  1) 
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possessed  the  symmetry  is  aad  V£  as  was  ©rigiaally 
assumed.  Satrisg  this  coaformattieo  of  toe  assumptions  made 
os  toe  field,  the  method  was  cased  to  geaerate  the  ieitial 
conditions  required  by  the  method  of  lines  solution  using 
5,  4,  5,  6,  zsd  S  lines.  For  the  subcritical  airfoil  the 
resultisg  Cp  di  stribut  j"--s  are  compared  to  each  other  is 
Fig.  7  on  the  next  page.  The  distribatios  obtained  from 
the  S-line  solac.oc  was  slightly  lower  than  that  obtained 
with  the  6- line  solution  but  they  were  so  close  that  the 
experimental  d at;,  could  not  be  read  accurately  enough  to 
determine  which  '-'as  the  better  solution.  Consequently,  it 
was  concluded  that  6- lines  were  sufficient. 

The  conclusion  which  can  be  drawn  from  these  results  is 
that  the  method  of  lines,  given  the  proper  initial  condition, 
produces  very  accurate  solutions  with  low  order  approximations. 
Further,  since  only  12.2  seconds  of  CDC  6600  computer  time 
was  required  tc  generate  all  the  Uj ,  Vj ,  velocity  ratios, 

and  CD  distribution  for  the  6- line  case,  the  solution  techni¬ 
que  satisfies  the  time  restraint  objective  required  of  the 
uethod . 

Critical  and  Supercritical  Airfoils 

_ _ — .  ■■  —  -  ■  -  *-  -  -  —  -  -  -- 

Kith  the  approach  established  on  the  subcritical  airfoil 
the  technique  was  next  tested  at  the  critical  Mach  number 
for  this  airfoil,  =  0.832.  Using  the  6-line  approximation 
the  results  were  found  to  be  in  excellent  agreement  with 
the  experimental  data;  the  C  distribution  was  again 
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coincidental  i to  tfee  ezp«?xinbeaitaB  values  and  at  tfee  grid! 
posau  z  *  ©.5©2.  y  =  O.fflSS*  cie  3C2cn  tunafeer  was  found  to  fee 
0.99S2  a5  ccctparecS  t©  tfee  ezperisneatalBy  predicted  vaBee  off 
M;  =  1  at  the  point  z  =  ©_5,  y  =  ©.©3. 

For  a  tfeird  rest  of  tfee  solution  tecfenique  tfee  super¬ 
critical  case  at  =  ©.S3  was  seEected  becaose  it  sr£BI 
appears  to  meet  tfee  symmetry  required  in  predicting  tfee 
initial  conditions-  As  discussed  in  Section  13,  tfeere  is  a 
cfeaage  of  signs  in  the  e s-Q  rel2ti©nsfeip  given  in  Eq  (22) 
wfeenever  a  sonic  line  is  crossed  in  tfee  flow  field,  Tfee 
computer  program  was  able  to  determine  tfee  points  at  wfeicfe 
tfeis  cfeange  W2S  to  occ^r  and  properly  predicted  tfee  sign  out 
tfee  reflated  requirement  for  tfee  slope  of  2^  to  fee  zero  at 
tfee  same  points  (refer  to  Figs.  4  and  5)  introduced  22 
additional  complication.  Consider  tfee  differential  equation 
for  0£ 


30J  f*(l  -  i/!I) L  SSj  (1  -  i/!s) 2L  3ri 

3£  (L  *  f)(l  -  jC|)2  3 H  (L  *  f)(l  -  |S|)2  3«1 

Letting  £*  be  the  value  at  which  the  sign  change  is  to  occur, 

then  over  a  snail  region  surrounding  this  point,  £  =  £*  ±  e 

3vi 

the  values  of  coefficients  and  nay  be  considered  constant 

and  this  expression  is  siaplified  to 

32  -  32- 

?T"  =  ci  sr  -  c2 
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Therefore,  for  a  change  am  sign  t©  ©ccer  im  ■J-^—  arithim 

*  S><j, 

tribe  region  £*  ♦  £  b  tr  can  fee  seen  that  a  very  accurate 

ap  ^roximation  t©  is  recoSred.  Trying  6  Biases  in  trie 

solation  techmicee  ([which  csed  2  second  order  approximation 

for  tribe  deriwatiTe}  failed  tr©  fee  accurate  emocagh  t©  get  a 

3f£= 

change  ias  sign  ©f  ~  when  the  sign  change  ©ccorred  im  the 

SC 

©-S  relationship  which  implies  that  2  larger  member  of  lines 
is  needed-  This  wocld  not  pose  any  serioas  complication 
except  that  tribe  met  feed  csed  t©  determine  initial  conditions, 
although  efficient  for  2  small  member  of  lines.  Deceases  2 
rerr  tedioas  2nd  lengthy  process  if  more  terna  eight  lines 
ore  rcqoired.  Since  there  was  no  war  to  determine  beforehand 
how  m2ay  lines  won  Id  be  repaired  to  get  the  accuracy  needed 
2nd  because  of  the  time  constraint  on  the  study,  the 
derelopmeat  of  this  solution  had  to  be  set  aside  in  order 
to  worh  on  more  generalized  problems. 


Initial  Conditions 

By  considering  only  problems  vith  symmetrical  airfoils 
at  zero  incidence  angle  and  by  restricting  ii  to  values  of 
i)  >  0,  it  will  be  recalled  that  the  boundary  conditions 
corresponding  to  the  infinity  conditions  are 


uj  :  vj  =  !!j  =  0  if  ^  =  ♦  1  or  n  =  1 


(27) 


while  at  the  boundary 


f'  if  0  <  (  < 


v0  = 


1  +  L 


(28) 


0  if£<0or£> 


1  ♦  L 
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Slmce  cmt-  soBtst  I©m  tietrBan sates  cases  a  auer-.-jiaas  roetlSxsflB  am  c&e 
£— direct bcce  TTd  hc^ax£xzy  comditlems  sas  available  me  botm 
£  =  i  B„  cMs  scgtricsts  7&jmt  the  ImiuiaB  cesditioms  are  wtBS 
established  the  solatiom  cam  progress  1ms©  difficulty 
becaese  ©f  two  comp II cat I ©ms  which  resoBt  freta  reBatlmg  the 
booaidary  xm£  Imitial  ccmditioss  to  the  real  physical  probBeaj. 
Flrst,  If  the  solce  a  ©a  was  started  from  5  =  -1  amd  the 
imitial  comditloss  were  aSB  set  eana B  t©  aero  as  predicted  by 
the  feoamdary  coadltioES  zhe  see  of  ©rcisEry  differential 
eooatioES  weal©  all  canal  zero  aad  Zhe  techaiaae  womld 
predict  aa  csdi sterbe©  flow  everywhere  cpseream  of  zhe  leading 
edge,  Slace  It  Is  knows  that  the  presence  of  toe  airfc-Il 
does  disturb  zhe  flow  upstream,  If  toe  Telocity  Is  less  teas 
supersonic.  It  Is  accessary  to  start  toe  solution  from  a 
point  away  frees  tie  boundary  £  «  -I  where  the  Initial  condi¬ 
tions  are  not  all  zero.  Second,  since  In  general  f*  ~  0  at 
the  trailing  and  leading  edge  of  an  airfoil,  the  boundary 
condition  at  ii  =  0  causes  discontinuities  in  vq  at  each  of 
these  points.  If  it  is  attempted  to  natch  past  these  points 
the  jusp  in  value  of  Vq  is  propagated  throughout  the  set  of 
governing  equations  and  the  system  is  found  to  be  unstable. 

To  avoid  this  instability  attempts  were  made  to  slightly 
alter  the  shape  of  the  airfoil  at  the  leading  and  trailing 
edges  such  that  Vq  changed  fron  0  to  f*  over  a  snail  region 
rather  than  at  a  single  point. 

In  developing  this  theory  for  predicting  the  initial 
condition  snoothing  functions  of  the  type 
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Tg  =  f«{£j  |  ;  0  <  s  <  £ 

V0  =  feS©HB  *  |);  -£  <  2  <  © 
t0  =  f - (£)  |  *  |c  :  '£  -  z  -  £ 

tq  =  A  *  c©s  *(1  +  |)1;  0  5  X  £  £ 

were  zest.eS  at  the  leadisg  edge  while  the  initial  coaditioss 
at  £  =  -0.9  were  perturbed  is  z  rariety  of  ccmbiaatioas  zbS 
magsitedes.  Sfith  this  eats  it  «ss  hoped  to  correlate  the 
change  in  initial  conditions  to  the  change  in  the  solution 
©rer  the  airfoil  bat  no  combination  of  initial  conditions 
and  smoothing  functions  could  be  found  to  give  stable  solu¬ 
tions  past  the  leading  edge.  The  solution  over  the  sub- 
critical  airfoil  was  produced  by  the  method  described 
previously  at  this  time  and  with  this  4-line  solution  in 
conjunction  with  the  forcing  function 

v o  =  f * (0.9) [(1.0  -  x)/0. 1] 

a  concentrated  effort  was  nade  to  narch  rt'f  the  back  of  the 
airfoil  and  correlate  the  change  in  initial  conditions  at 
x  =  1/2  to  the  flow  field  that  resulted  at  x  =  2.  It  was 
found  that  a  perturbation  on  the  order  of  10'^  inposed  on 
any  one  of  the  initial  conditions  would  substantially  alter 
the  flow  at  x  =  2  and  that  although  the  nagnitudc  of  all  the 
u £  and  decreased  downstream  of  the  trailing  edge  they 
eventually  would  diverge  and  no  set  of  initial  conditions 
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could  be  fooisd  co  force  convergence  to  undisturbed  f lev  far 
downstream.  Because  of  the  sensitivity  of  tbe  flow  field 
to  tbe  initial  conditions,  tbe  value  of  a?  and  v^  marching 
oa  or  off  both  ends  of  tbe  smoothing  functions  are  critical 
to  tbe  solution.  Consequently,  no  way  could  be  found  to 
separate  tbe  influence  of  a  change  in  initial  conditions 
froas  tbe  error  injected  by  tbe  smoothing  functions  and  with¬ 
out  this  separation  no  rationale  could  be  found  for  the 
selection  of  either  one. 

Ifavy  5Lfall  Solutions 

It  was  realized  very  late  in  the  study  that  the  solution 
technique  could  be  applied  equally  well  to  an  infinite  wavy 
wall.  The  exact  solutions  for  Cp,  u,  and  v  for  both  the 
subsonic  and  supersonic  case  are  developed  by  Liepnann  and 
Roshko  (Ref  6}  and  offer  a  very  efficient  neans  of  testing 
the  proposed  solution  technique.  Although  these  exact  solu¬ 
tions  cannot  be  used  to  deternine  the  applicability  of  the 
nethod  to  transonic  problems,  they  can  serve  to  evaluate  the 
accuracy  of  the  solution  technique  in  predicting  the  flow 
field  as  given  by  the  exact  solutions  if  the  method  is 
started  with  the  correct  initial  conditions. 

To  perform  this  evaluation  the  system  of  equations  were 
first  reduced  to  the  simplified  problem  solved  by 
Liepmann  and  Roshko  by  approximating  the  ft-u  relationship 
given  in  Eq  (22)  as 

ft  =  (1  -  M  2)u 

v  GO  * 
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and  by  simplifying  the  y-n  coordinate  transforaation  to 


n  = 


y  *  l 


so  the  boundary  condition  at  y  =  f(x)  could  be  approxinately 
applied  at  y  =  0.  The  wall  function,  f (x) ,  was  defined  by 

f(x)  =  0.01  sin  2ttx 

where  the  wave  aaplitude  was  selected  as  0.01  with  a  wave 
length  of  1.  By  choosing  x  =  0  as  the  starting  point  the 
initial  conditions  were  defined  as 


ui 


=  0 


vi 


=  (0.01)  (2it) e 


rtt  (2»>(l-Mw)1/2 


i  =  0,  1,  2,  ...,  (N-l) 


for  the  subsonic  case  and 


Ui  =  -  (0.01)(2n)(Mco2-l)"1/2  cos  [-2ir  M^2-!)172] 


Vi  =  (0.01)  (2tt)  cos  [ - 2 it  (Mro2-1)  1/2] 

N-i 


i  =  0,  1,  2,  ....  (N-l) 


for  the  supersonic  case  where  N  is  the  number  of  lines  to 
be  used  in  the  solution.  Finally,  in  order  to  test  the 
technique  on  both  cases,  Mro  =  0.2  and  M  =  2.0  were  arbitrar¬ 
ily  selected  as  the  subsonic  and  supersonic  Mach  numbers. 
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Because  finite-differences  are  used  to  approximate  the 
derivatives  and  since  u  and  v  are  changing  most  rapidly  in 
the  vicinity  of  the  wall,  it  is  necessary  to  have  a  good 
distribution  of  lines  in  this  region  if  the  solution  is  to 
be  accurate.  Solving  Eq  (29)  for  y  and  letting  r\  -  i/N 
implies  that 

iL 

^i  “  N  -  i 

from  which  it  can  be  observed  that  the  value  of  L  determines 
how  close  the  ith  line  is  to  the  body.  With  a  large  number 
of  lines  (N  >  15)  many  of  the  lines  are  naturally  close  to 
the  body  due  to  the  coordinate  transformation  and  a  change 
in  L  (0  <  L  <  3)  has  little  effect  on  the  solution.  For  a 
smaller  number  of  lines  I.  becomes  increasingly  important  and 
for  less  than  8  lines  a  substantial  change  in  the  solution 
was  observed  as  L  was  varied  between  0  and  2.  Figure  8  on 
the  page  following  this  section  shows  the  4-line  solution 
for  the  subsonic  wavy  wall  with  L  »  0.5  while  Fig.  9  shows 
the  improvement  in  this  solution  as  L  is  changed  to  0.65. 

In  general,  an  optimum  value  of  L  could  be  found  for  every 
value  of  N  but  in  each  case  tested  L  =  0.5  gave  reasonable 
results.  Therefore,  L  =  0.5  was  selected  as  the  reference 
value  of  l  used  in  all  solutions. 

As  was  pointed  out  in  the  formulation  of  this  solution 
technique,  the  boundary  condition  at  infinite  was  defined  to 
be 

u^  =  v^«=fij,=0  0  |£|  =  1  or  n  =  1 
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x/c 


Fig.  8.  Solution  Obtained  for  the  Cp  Distribution  over  a 
Subsonic  Wavy  Wall  Using  4  Lines  witn  L  =  0.5  as  Compared 
to  the  Exact  Solution. 
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For  the  subsonic  wavy  wall  the  exact  solution  agreed  with  this 
condition  and  consequently  the  solution  required  only  4-lines 
to  get  the  accuracy  as  shown  in  Fig.  9.  However,  with  the 
supersonic  wavy  wall  the  exact  solutions  predicts  that  u  and 
v  vary  sinusoidally  at  infinity  and  although  this  solution 
gives  the  range  of  these  values  it  can  not  be  used  to  predict 
specific  values  for  u  or  v  at  particular  points  on  the 
boundary.  When  u,  v,  and  A  were  set  to  zero  at  n  *  1  errors 
resulted  in  the  solution.  However,  as  the  number  of  lines 
was  increased  the  solution  converged  rapidly  and  in  Fig.  10 
the  Cp  distribution  obtained  using  20  lines  is  shown  to 
compare  very  well  with  the  exact  solution. 


32 


GAM/AE/72-7 


IV.  Conclusions  an  &  Pcccatasesdiations 

Conclusions 

From  the  results  obtained  for  the  subsonic  and  super¬ 
sonic  wavy  vail  and  the  solutions  found  for  the  circular  arc 
airfoil  it  can  be  concluded  that  the  proposed  solution 
technique  is  capable  of  producing  extrenely  accurate  results 
with  very  short  computer  tiaes  provided  it  is  furnished  with 
correct  initial  conditions. 

Recoanendations 

The  following  recoanendations  are  proposed  for  future 

work : 

1.  It  is  suggested  that  for  this  problen  an  unsteady 
systen  of  governing  equations  should  be  used.  With  this 
three-dimensional  coordinate  system  the  marching  method  can 
be  set  up  in  the  time  direction  and  a  method  of  weighted 
residuals  could  be  used  to  approximate  the  derivatives  in 
the  x  and  y  directions.  In  this  form,  the  problem  of 
determining  generalized  initial  conditions  can  be  separated 
from  the  problem  of  handling  the  discontinuities  in  the  x-y 
plane  and  their  individual  effects  can  be  evaluated  and 
controlled . 

2.  A  coordinate  transformation  of  the  form  used  in  this 
development  should  be  used  to  control  the  domain  of  the 
independent  variables.  Although  the  opportunity  was  not 
available  to  demonstrate  the  advantage  of  shrinking  the 
airfoil  thickness  to  zero,  it  is  suggested  that  this  benefit 
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be  retailed  because  of  the  thick  airfoils  used  is  the 
transonic  regime. 

5.  The  introduction  of  2  to  separate  toe  nonlinearity 
from  tee  differential  equation  made  the  system  much  easier 
to  sfork  vith  and  it  is  suggested  that  this  same  technique 
should  be  applied  whenever  similar  nonlinear  terms  appear 
in  the  differential  relations. 

4.  The  method  of  lines,  as  used  in  this  development, 
is  highly  recommended  as  a  technique  for  reducing  the  partial 
differential  equations.  Besides  being  a  *Tery  straightforward 
nethod  to  apply  it  has  an  outstanding  advantage  in  that  the 
resulting  unknowns  can  be  easily  identified  vith  the  physical 
problen. 
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Appendix  A 

EtetaA  5  s  of  tbe  So  Bos  ice  Tecbaicgoe  Forme Saties 

This  appendix  serves  as  as  exteasios  to  the  develepmeat 
as  preseated  is  Sectioa  II  of  the  maia  bodjr.  Its  parpose 
is  to  give  a  detailed  discession  and  formalatioa  of  the 
important  mathematical  steps  in  areas  where  details  of  tbe 
solution  technique  might  be  desired  by  some  readers.  Three 
primary  subjects  will  be  covered:  (1)  tbe  aondimensionaii- 
zation  of  the  basic  governing  equations,  (2)  tbe  coordinate 
transformation,  and  (5)  tbe  application  of  tbe  method  of 
lines. 


Kondimensionalization 

Tbe  set  of  governing  equations  which  were  selected  for 
this  presentation  are  the  small  perturbation  relationship 


u  -  h^)  ,  |X 


3x 


oy 


2  y+l 


u 


3u 

3x 


which  nay  fee  written  by  rearranging  the  derivatives  as 


3 

Jx 


1(1  -  He 


.  - 


H. 


2  Hi 

»  U„ 


3v 


-}+  jy 


=  o 


(29) 


plus  the  irrotationality  condition 


3v  3u 
Tx  "  3~y 


0 


(30) 


In  these  equations  and  arc  the  free  stream  Mach  number 
and  uniform  velocity,  y  is  the  ratio  of  specific  heats,  and 
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o  aad  r  are  the  tpertorbatioo $  vhi ch  are  def£ffles3  by 

0  =  l®„  *  a}es  *  re? 

where  IF  is  the  local  total  Telocity  rector.  These  differen¬ 
tial  goremiag  eqeatioas  are  sabjiect  to  the  boasdary  cosdi- 
tioES  at  x  or  y  eacals  iaficity  where 

o  =  r  =  0  (31) 


aad  on  tbe  serface  of  the  body  where 


r  _  df 
a  ♦  !!,_>  "  dx 


(32) 


2nd  y  =  f(x)  defines  the  shape  of  the  airfoil- 

To  generalize  these  relationships  it  is  desirable  to 
express  therm  in  their  nondimensional  form.  To  do  so  the 
lengths  are  nondimensionaii zed  with  respect  to  the  chord 
length 


x*  = 


x 

c 


y 

c 


f*  = 


f 

c 


and  the  velocities  are  nondinensi onalized  with  respect  to 
the  free  stream  velocity 


u*  = 


u_ 

U. 


v* 


V 

U«> 
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SOaiisg  zhese  defiasitioas  aadi  safes t£ to ting  into  Eqs  (29) 
tferoa gsi  (32)  cfeaeges  toe  form  of  these  espressiosis  to 


1(1  -  S«as2)a*  -  4  M^CT+ljc*2!  ♦  =  0  C33) 


3v* 

5s* 


3s1 

Jp 


=  0 


(34) 


a*  =  v*  =  0  §  <= 


(35) 


df 


r* 


a*  *  1 


di1 


£  v*  = 


(36) 


Dropping  the  superscripts  for  simplicity  of  notation  gives 
the  equations  2S  the}*  are  used  in  the  main  body. 

Coordinate  Trans  formation 

It  was  decided  to  use  a  coordinate  transformation 
because  it  is  desirable  to  limit  the  domain  of  the  indepen¬ 
dent  variables.  The  specific  transformation  chosen  to 
accomplish  this  objective  is; 


s  = 


M  -  L 


(37) 


n  =  y  '  f(5° 


|y|  -  l 


(38) 


where  L  is  some  positive  constant.  The  relationship  between 
the  original  independent  variables  and  the  new  coordinates  is 
shown  graphically  in  Figures  A1  and  A2.  As  can  be  seen  from 
these  figures,  as  the  Jx|  or  |y|  go  to  infinity  the  |£|  and 
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Fig.  A2.  Relationship 
between  the  Independent 
Variables  v  -  q. 


Fig.  Al.  Relationship 
between  the  Independent 
Variables  x  -  £- 

fill  are  restricted  to  values  of  less  than  one  therefore 
these  transformations  reduce  the  douain  of  the  independent 
variables  to  finite  regions.  Also  Fig.  A2  shows  the  addi¬ 
tional  benefit  of  reducing  the  thickness  of  any  airfoil  in 
the  x-y  plane  to  zero  in  the  £-n  plane  but  it  should  be 
noted  that  the  plot  of  y-n  is  dependent  on  the  value  of  x 
because  f(x)  appears  in  Eq  (38). 

In  order  to  transforn  the  governing  equations  into 
the  new  coordinate  system  the  partial  derivatives  must  be 
expressed  in  terns  of  £  and  n .  If  some  dependent  variable 
is  defined  by 


g(x,y)  =  <U£,n) 

then  the  partial  derivatives  with  respect  to  x  and  y  may  be 
written  by  the  chain  rule  as 
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3g  _  3G  3£  3G  3  si 
3x  “  ®C  ex’  +  35)  «x 


(39) 


pg  3G  SC  3G  35) 

3y  -  Sc  ey  +  Sn  37  (40) 

Fro®  Eos  (37)  and  (38)  the  partials  of  C  and  n  vith  respect 
to  x  and  y  are  evaluated  to  be 


(M  ♦  l)2 


3£ 

3?=  ° 


35)  _  -f* 

31  '  II  I  *  L 


3n  =  l  ♦  1  f  ( 

(|yj  ♦  L)2 


or  in  terns  of  the  new  independent  variables 


H  (1  -  Ul)2 

3x  ’  l - 


=  -f<il  '  M) 

L  ♦  Ifl 


(1  -  ini)2 


L  ♦  f| 


For  flow  fields  symmetrical  about  n  =  0,  only  the  upper 
half  plane  is  considered  and  Hqs  (35)  and  (40)  may  be 
written  as 
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is  _  Cl  -  Ul)2  3G  f  *  (1-n)  3G 
3x  L  5T  L+f  7rf 


(41) 


3g  _  (l-n)2  3G 
Ty  ~  L+f  off 


(42) 


Returning  to  the  differential  governing  equations  as 
given  in  Eqs  (33)  and  (34)  and  defining 

A  =  (1  -  Ma92)u  -  1  M002(Y+1)u2  (43) 

then  with  the  definitions  given  in  Eqs  (41)  and  (42)  the 
governing  equations  in  the  new  coordinate  system  are 

(i  -  Ul)2  M  f " (l-n)  3«  +  (l-n)2  3v  .  0 

L  9  5  L+f  Sn  L+f  3n 

(i  -  Ul)2  3v  f  *  (l-n)  3v  (l-n)2  3u  _ 
l  35  L+f  an  "  L-f  5n"  "  u 


To  get  the  form  of  these  equations  desired  for  application 
of  the  method  of  lines  solution  technique,  derivatives  with 
respect  to  5  remain  on  the  left  hand  side  such  that 


9fl 

f  •  (l-n)L 

BQ  _ 

(l -n) 2l 

9v 

3C  ' 

(L+f)  ( 1  -  Ul)2 

3n 

(L+f)  (1  - 

Ul)2  3T1 

9v 

f ’ (l-n) L 

+ 

(l-n) 2l 

9u 

H  = 

(L+f)  (1  -  |5|)2 

3n 

(L+f)  (1  - 

Id)2  Sn 

(44) 


(45) 


Finally,  the  boundary  conditions  as  given  in  Eqs  (35)  and 
(36)  become 

u  =  v  =  fi  =  0  0 | 5 | or  n  =  1 
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and 


=  f'  e  P  =  0 

u+ 1 

< 

in  the  new  coordinate  system. 

App lying  the  Method  of  Lines 

In  general,  the  functions  u  and  v  are  changing  much 
more  rapidly  in  the  ^-direction  than  in  the  p-direction. 
Therefore,  to  use  the  method  of  lines  the  derivatives  with 
respect  to  n  are  chosen  to  be  represented  by  a  finite- 
difference  approximation.  Along  any  line  parallel  to  the 
p-axis,  except  for  p  =  0,  the  first-order  partial  of  a 
function  G(£,p)  with  respect  to  p  is  approximated  by 

3G  „  G($,p  4  Ap)  -  G  (£,  p  -  Ap) 

sp - mr -  (46) 

which  is  the  second-order  central-difference.  For  the  line 
p  =  0,  since  only  values  of  p  greater  than  zero  are  con¬ 
sidered,  the  partial  is  approximated  by  the  second-order 
forward  difference 


3G 

ap 


p  =  0 


~  t-3G(5,0)  4  4G (£ , Ap)  -  G(C,2An)] 


(47) 


where  in  both  Hq  (46)  and  (47)  the  Ap  refers  to  the  distance 
between  evenly  spaced  lines.  Dividng  the  p-domain  by  (N-l) 
of  these  parallel  lines  between  p  =  0  and  p  =  1, 


P  =  0, 


i  1 

N  *  N  ’ 
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then  the  value  of  n  along  the  ith  line  is  written 


(48) 


Introducing  subscripts  corresponding  to  the  value  of  n  and 
using  Eq  (48)  implies  that  the  partials  in  Eqs  (46)  and  (47) 
may  be  written 


1£  _  N  , 

3q  ~  2  (Gi+1  '  Gi-lJ 


(49) 


3G0  ^ 

=  2  f’3G0  +  4G1  "  G2^  (50) 

Applying  these  definitions  in  the  governing  equations  as 
given  in  Eqs  (44)  and  (45)  reduces  these  partial  differential 
relationships  to  a  set  of  first-order  ordinary  differential 
equations  of  the  form 


dni  =  f'(l  -  i/N)L  M 
dC  (L+f)  (  1  -  Ul)2  2 


Vl 


)  - 


(1  -  i/N)2L  N 
(L+f)  (1  -  |S  |)  2  2 


(vi  +  l 


vi-l) 


dflo  =  QL _ ,  o 

dS  2  (L+f)  Cl  -  U  j  )  2  ° 


4ftj  •  fl2)  - 


LN 


2 (L+f ) (1  -  |S|) 


(-3v, 


4v. 


V2} 
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Appendix  B 


A  Solution  Technique  Based  on  the  Method 
of  Weighted  Residuals 


This  appendix  presents  an  alternate  solution  technique 
for  the  transonic  airfoil  problem.  This  is  the  original 
approach  developed  for  this  study  but  had  to  be  abandoned 
when  no  rationale  procedure  for  initializing  the  variables 
could  be  found.  Since  the  solution  technique  appears  power¬ 
ful  enough  to  warrant  further  investigation  if  the  initiali¬ 
zation  problem  can  be  overcome,  a  complete  formulation  of 
the  theory  is  given  here. 

The  set  of  governing  equations  used  by  this  technique 
are  the  same  as  those  used  with  the  method  of  lines.  After 
nondimensionalizing  and  applying  the  coordinate  transforma¬ 
tion  the  set  consists  of 


dQ 

n 


Lf 1  (i-n)  da  L(i-n) 2  9_v 

2  **  (L+fj ( i  -  u.!r  '3n 


(L+f)d  -  id)' 


(51) 


9  v 


Lf(l-n)  9  v 


L(i-n) 


9u 


(L*f)(l  -  Uj)2  9n  (L+  f )  ( 1  -  U|)2  9n 


(52) 


u  = 


a  -  ) 

Mco2(Y  +  1) 


1  ♦  /  1  - 


2M„qY+i) 


2>2 


a 


q-M» ) 


(53) 
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The  boundary  conditions  to  be  satisfied  are:  at  the  body 
surface,  rj  =  0, 


Sv  _  _L _  f„ 

35  (1  -  Ui)2 

and  at  infinity,  £  or  n  =  +1, 

■  u  =  v  =  ft  =  0 


,  (54) 


(55) 


For  this  formulation  the  dependent  variables  u,  v,  and  ft  are 
approximated  by  a  series  expansion  in  n  where  the  terms  are 
arranged  such  that  the  series  automatically  satisfy  the 
boundary  condition  at  n  =  1 

M 

u(£,n)  =  l  ai(?)(l-n)1  (56) 

i=  1 

N  , 

v(f,n)  =  l  bi(C)d-n)J  (57) 

3*1' 

'  p  t 

«(5.n)  =  I  ck(C)  ( l -n) ■  '  (58) 

k=l 

Substituting  these  into  Eqs  (49)  through  (52) 


l  ck' (l-n) 


-L 


k=  1 


(L  +  f)(l  -  U|) 


f*  l  ckk(i-n)k  -  l  bU(i-n)j  +  1 


k=  1 


3  =  1  ' 


(59) 
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N 

l  b-'d-n)3 
j=i  3 


_ -L _ 

(L  ♦  f)(l  -  is!)2 


I  bJci-n)j  +  ?  aj1  (l-n) i+1 

j=l  i=l 


(60) 


M 

X 


i  =  l 


aiCl-n)1  = 


(1  ~  »J) 

M«o2  (Y  +  1) 


1  ± 


2mJ(Y  +  1)  P 
(1  -  M J)2  k=l 


ckd-n)k 


(61) 


N 


I 

j  =  l 


_ L 

(1  -  Ul)2 


f" 


(62) 


In  this  set  of  equations  at  a  given  n  there  are  a  total  of 
(M  +  N  +  P)  unknowns;  namely,  all  the  coefficients  ai ,  bj, 
and  ck.  Therefore,  a  unique  solution  for  this  sytem  requires 
that  the  same  number  of  independent  equations  be  derived 
from  Eqs  (59)  through  (62).  The  method  of  weighted  residuals 
offers  many  techniques  by  which  this  system  of  equations  may 
be  developed  (Ref  3).  Of  these  techniques,  the  collocation 
method  is  the  easiest  to  apply  and  was  selected  for  use  here. 
Basically,  this  method  assumes  that  an  approximate  solution 
for  the  coefficients  a^,  b^,  and  Cj.  can  be  obtained  if  the  Eqs 
(59)  through  (62)  are  satisfied  at  a  finite  number  of  n-values. 
The  collocation  method  allows  the  user  to  select  the  number 
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of  each  type  of  equation  and  the  particular  values  of  n  to 
be  used  in  obtaining  these  equations.  In  this  application 
the  equations  were  selected  in  the  following  Banner.  By 
allowing  n  to  take  on  the  values 


in  Eq  (59),  P  ordinary  differential  equations  are  obtained 
involving  the  cj.1  terns.  In  Eq  (60),  n  is  defined  to  be 

0  =  |  ;  „  =  1,  2,  ...,  (N-l) 

which  together  with  the  boundary  condition  in  Eq  (62)  gives 
N  ordinary  differential  equations  involving  the  b j '  terns. 
Finally,  M  algebraic  relationships  are  obtained  from  Eq  (61) 
by  letting  n  take  on  the  values 

»  -  *  1,  2 . M 

Collecting  this  set  gives  the  complete  new  system  of  governing 
equations 
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(66) 


The  sign  choice  in  Eq  (66)  is  made  in  the  same  manner  as 
discussed  in  Section  II  of  the  main  body  where  the  critical 
value  of  u  for  this  case  is 


To  solve  this  system  observe  that  if  initial  conditions 
are  known  for  the  coefficients  bj  and  Cj.  then  Eq  (66)  can  be 
solved  for  the  corresponding  values  of  a^  .  Everything  is 
then  known  on  the  right  sides  of  P.qs  (63)  through  (65)  and 
the  system  of  ordinary  differential  equations  can  be  solved 
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by  the  Kuage-Sutta  method  to  set  up  the  marching  method  used 
throughout  the  field.  Unfortunately,  no  rationale  coaid  be 
found  for  determining  the  initial  conditions  needed  to  start 
this  solution  technique.  Sfitn  the  method  of  lines  the 
dependent  variables  to  be  determined  were  the  perturbations 
at  various  values  of  q  (u^  and  v^) .  A  knowledge  of  typical 
distributions  of  u  and  v  then  could  be  used  to  determine 
the  required  initial  condition  when  the  solution  was  started 
on  the  airfoil  surface.  Kith  this  collocation  technique 
though,  the  variables  to  be  solved  for  are  the  coefficients 
of  series  (a^,  bj,  and  cj-)  and  since  it  was  unknown  how  a 
typical  distribution  of  these  terns  night  look,  it  was  not 
possible  to  use  a  sinilar  rationale  for  determining  the 
initial  conditions. 


51 


GAM/AE/72-7 


Vita 

John  Richard  McCracken  was  born  on  11  March  1944  in 
Denver,  Colorsdo.  He  entered  the  Air  Force  on  13  March  1961, 
attending  the  University  of  Maryland  and  then  Auburn 
University  from  which  he  received  the  degree  of  Bachelor  of 
Aerospace  Engineering  in  March  1970  and  was  elected  to  Tau 
Beta  Pi  and  Signa  Ganna  Tau.  After  receiving  his  commission 
in  the  USAF  in  June  1970  he  was  assigned  to  the  Air  Force 
Institute  of  Technology. 

Pernanent  address:  Box  355 

Hayden,  Colorado  81639 


This  thesis  was  typed  by  Jane  Manemann. 


52 


